﻿using RayDen.Library.Core.Primitives;

namespace RayDen.Library.Core.Noise {
    public class NoiseProvider {
        public static int NoiseDim = 15;

        private static NoiseProvider instance;

        public static NoiseProvider Instance {
            get {
                if (instance == null)
                    instance = new NoiseProvider();
                return instance;
            }
        }

        private float[, ,] NoiseTable;


        NoiseProvider() {
            InitNoise();
        }

        public static float Spline(float t) {
            return t * t * (3f - 2f * t);
        }

        public static float SineWave(float x) {
            return 0.5f * (1.0f + MathLab.Sin(x));
        }

        public static float Mod(float x, float y) {
            if ((x = fmod(x, y)) < 0)
                return x + y;
            else
                return x;
        }

        private static float fmod(float x, float y) {
            return x % y;
        }

        public void InitNoise() {
            var RNG = new FastRandom();
            NoiseTable = new float[NoiseDim + 1, NoiseDim + 1, NoiseDim + 1];
            int i, j, k;
            for (i = 0; i < NoiseDim; i++)
                for (j = 0; j < NoiseDim; j++)
                    for (k = 0; k < NoiseDim; k++)
                        NoiseTable[i, j, k] = RNG.NextFloat();
        }

        public float Noise(Vector p) {

            float sx = Mod(p.x, NoiseDim);
            float sy = Mod(p.y, NoiseDim);
            float sz = Mod(p.z, NoiseDim);
            int ix = (int)sx;
            int iy = (int)sy;
            int iz = (int)sz;
            int jx = ix + 1;
            int jy = iy + 1;
            int jz = iz + 1;

            if (jx >= NoiseDim)
                jx = 0;

            if (jy >= NoiseDim)
                jy = 0;

            if (jz >= NoiseDim)
                jz = 0;

            sx = Spline(sx - ix);
            sy = Spline(sy - iy);
            sz = Spline(sz - iz);

            return (1 - sx) * (1 - sy) * (1 - sz) * NoiseTable[ix, iy, iz] +
                   (1 - sx) * (1 - sy) * sz * NoiseTable[ix, iy, jz] +
                   (1 - sx) * sy * (1 - sz) * NoiseTable[ix, jy, iz] +
                   (1 - sx) * sy * sz * NoiseTable[ix, jy, jz] +
                   sx * (1 - sy) * (1 - sz) * NoiseTable[jx, iy, iz] +
                   sx * (1 - sy) * sz * NoiseTable[jx, iy, jz] +
                   sx * sy * (1 - sz) * NoiseTable[jx, jy, iz] +
                   sx * sy * sz * NoiseTable[jx, jy, jz];
        }

        public float Noise(float x, float y, float z)
        {

            float sx = Mod(x, NoiseDim);
            float sy = Mod(y, NoiseDim);
            float sz = Mod(z, NoiseDim);
            int ix = (int)sx;
            int iy = (int)sy;
            int iz = (int)sz;
            int jx = ix + 1;
            int jy = iy + 1;
            int jz = iz + 1;

            if (jx >= NoiseDim)
                jx = 0;

            if (jy >= NoiseDim)
                jy = 0;

            if (jz >= NoiseDim)
                jz = 0;

            sx = Spline(sx - ix);
            sy = Spline(sy - iy);
            sz = Spline(sz - iz);

            return (1 - sx) * (1 - sy) * (1 - sz) * NoiseTable[ix, iy, iz] +
                   (1 - sx) * (1 - sy) * sz * NoiseTable[ix, iy, jz] +
                   (1 - sx) * sy * (1 - sz) * NoiseTable[ix, jy, iz] +
                   (1 - sx) * sy * sz * NoiseTable[ix, jy, jz] +
                   sx * (1 - sy) * (1 - sz) * NoiseTable[jx, iy, iz] +
                   sx * (1 - sy) * sz * NoiseTable[jx, iy, jz] +
                   sx * sy * (1 - sz) * NoiseTable[jx, jy, iz] +
                   sx * sy * sz * NoiseTable[jx, jy, jz];
        }

        public Vector Noise3D(Vector p) {
            Vector res = new Vector();

            float sx = Mod(p.x, NoiseDim);
            float sy = Mod(p.y, NoiseDim);
            float sz = Mod(p.z, NoiseDim);
            int ix = (int)sx;
            int iy = (int)sy;
            int iz = (int)sz;
            int jx, jy, jz;

            sx = Spline(sx - ix);
            sy = Spline(sy - iy);
            sz = Spline(sz - iz);

            for (int i = 0; i < 3; i++) {
                ix = (ix + 5) % NoiseDim;
                iy = (iy + 5) % NoiseDim;
                iz = (iz + 5) % NoiseDim;

                if ((jx = ix + 1) >= NoiseDim)
                    jx = 0;

                if ((jy = iy + 1) >= NoiseDim)
                    jy = 0;

                if ((jz = iz + 1) >= NoiseDim)
                    jz = 0;

                res[i] = (1 - sx) * (1 - sy) * (1 - sz) * NoiseTable[ix, iy, iz] +
              (1 - sx) * (1 - sy) * sz * NoiseTable[ix, iy, jz] +
              (1 - sx) * sy * (1 - sz) * NoiseTable[ix, jy, iz] +
              (1 - sx) * sy * sz * NoiseTable[ix, jy, jz] +
              sx * (1 - sy) * (1 - sz) * NoiseTable[jx, iy, iz] +
              sx * (1 - sy) * sz * NoiseTable[jx, iy, jz] +
              sx * sy * (1 - sz) * NoiseTable[jx, jy, iz] +
              sx * sy * sz * NoiseTable[jx, jy, jz];
            }

            return res;

        }


        public Vector Noise3D(float x, float y, float z)
        {
            Vector res = new Vector();

            float sx = Mod(x, NoiseDim);
            float sy = Mod(y, NoiseDim);
            float sz = Mod(z, NoiseDim);
            int ix = (int)sx;
            int iy = (int)sy;
            int iz = (int)sz;
            int jx, jy, jz;

            sx = Spline(sx - ix);
            sy = Spline(sy - iy);
            sz = Spline(sz - iz);

            for (int i = 0; i < 3; i++)
            {
                ix = (ix + 5) % NoiseDim;
                iy = (iy + 5) % NoiseDim;
                iz = (iz + 5) % NoiseDim;

                if ((jx = ix + 1) >= NoiseDim)
                    jx = 0;

                if ((jy = iy + 1) >= NoiseDim)
                    jy = 0;

                if ((jz = iz + 1) >= NoiseDim)
                    jz = 0;

                res[i] = (1 - sx) * (1 - sy) * (1 - sz) * NoiseTable[ix, iy, iz] +
              (1 - sx) * (1 - sy) * sz * NoiseTable[ix, iy, jz] +
              (1 - sx) * sy * (1 - sz) * NoiseTable[ix, jy, iz] +
              (1 - sx) * sy * sz * NoiseTable[ix, jy, jz] +
              sx * (1 - sy) * (1 - sz) * NoiseTable[jx, iy, iz] +
              sx * (1 - sy) * sz * NoiseTable[jx, iy, jz] +
              sx * sy * (1 - sz) * NoiseTable[jx, jy, iz] +
              sx * sy * sz * NoiseTable[jx, jy, jz];
            }

            return res;

        }
        public float Turbulence(Vector p, int Octaves) {
            float k = 1f;
            float res = 0;
            Vector r = p;

            for (int i = 0; i < Octaves; i++) {
                res += Noise(r) * k;
                r *= 2f;
                k *= 0.5f;
            }

            return res;

        }
        public Vector Turbulence3d(Vector p, int Octaves) {
            float k = 1f;
            Vector r = p;
            Vector res = new Vector(0f);

            for (int i = 0; i < Octaves; i++) {
                res += Noise3D(r) * k;
                r *= 2f;
                k *= 0.5f;
            }

            return res;

        }

        public Vector Turbulence(float x, float y, float z, int Octaves)
        {
            float k = 1f;
            Vector r = new Vector(x, y, z);
            Vector res = new Vector(0f);

            for (int i = 0; i < Octaves; i++)
            {
                res += Noise(r) * k;
                r *= 2f;
                k *= 0.5f;
            }

            return res;

        }
        public Vector Turbulence3d(float x, float y, float z, int Octaves)
        {
            float k = 1f;
            Vector r = new Vector(x,y,z);
            Vector res = new Vector(0f);

            for (int i = 0; i < Octaves; i++)
            {
                res += Noise3D(r) * k;
                r *= 2f;
                k *= 0.5f;
            }

            return res;

        }
    }

}
